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We study the physical mechanism of the surface d ^ d + is transition proposed as the interpre- 
tation of results of tunneling experiments into ah planes We base our argument on first-order 
perturbation theory and show that the zero-bias states drive the transition. We support the argu- 
ment by various estimates and consistency checks. 



I. INTRODUCTION 



It has now been firmly estabhshed that the order parameter A in homogeneous cuprate super- 
conductors has a (i-wave symmetry IQ. It follows that inhomogeneities can scatter quasiparticles 
between directions that experience opposite signs of A. This effect is strongest at a specularly re- 
flecting (110) surface, because A changes sign along each quasiclassical trajectory upon reflection 
(see Fig. 1). As a consequence of the Atiyah-Patodi-Singer index theorem 0, the Andreev equation 
along each such trajectory has in its spectrum a bound state at exactly zero energy, irrespective 
of the detailed shape of the potential 0. Collectively, these states then make up a peak in the 
tunneling spectra, which has been observed experimentally 
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FIG. 1. a) A schematic picture of the normal metal-superconductor junction in the (110) direction with a typical quasiclas- 
sical trajectory. 

b) A schematic graph of the pairing potential along the trajectory in a). 
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It was predicted theoretically [0] and confirmed experimentally [|] that below a certain temperature 
of the order of IK, this peak of the zero-bias states (ZBSs) splits. The theoretical interpretation 
of this splitting is that via a sub-dominant pairing interaction of a different symmetry, say s, a 
subdominant order parameter is induced close to the surface. It is phase-shifted by 7r/2 relative to 
the dominant (i-wave, which gives the total order parameter "c? + is" symmetry and indicates broken 
time-reversal symmetry. 

The full self-consistent calculations in the Eilenberger formalism that are based on this inter- 
pretation |]9|,p!0[| are in a good quantitative agreement with the experimental data. However, the 



mechanism of the transition is not manifest in the numerical solutions of the Eilenberger equations. 
We believe that understanding of the basic physics of the d d + is transition is especially needed 



now in light of recent experiments that call this interpretation into question |jTT]|. That is the purpose 
of the work presented in this paper. 

Below, we show by a simple argument based on first-order perturbation theory that the degrees of 
freedom driving the d d + is transition are the ZBSs, and that we can neglect the effect of all the 
remaining states. Hence, to understand the mechanism of the transition, we have to deal only with 
the ZBSs, which is convenient since these states are least sensitive to the unknown surface details. 

The paper is organized as follows: In Section |T|, we demonstrate our strategy on the familiar case 



of BCS instability. The main argument is presented in Section |TT| after we have extended the BCS 
formalism to inhomogeneous systems and non-s-wave pairing. Based on this argument, we calculate 
A at T = in Section |^ and estimate the transition temperature to the d + is state in Section 0. 



In Section VI, we discuss the surface current. Finally, we discuss our results in Section VII 



II. BCS INSTABILITY 

There are various ways to consider the energetic costs and benefits of the transition to the superfiuid 
state. The one that has proven useful in our study of the d ^ d + is transition is to decouple the 
attractive four-fermion interaction by the Hubbard-Stratonovich (HS) transformation, and to make a 
saddle-point (mean-field) approximation. That way, we break up the total free energy of the system 



2 



into free energy of single particle states, which is lowered by the gap A, and the extra term from the 
HS transformation, which grows (quadratically) with A. We then see that at small enough T, the 
system favors transition to the superfluid state. 
We will demonstrate this on the familiar BCS case. The model Hamiltonian is Q 



H = ^ek4,,Ck,a - \V\ 4|cLkiC^k'iCk'T 

k,(T 

which gives rise to the partition function 

f 



(2.1) 



/3 

■Jdr 



k,k' 



ko- k,k' 



(2.2) 



where the c's are now r— dependent Grassmann numbers. We perform the HS transformation by 
multiplying the partition function by the (infinite) constant 



so 



_ -\y\ J (i-r X] (<^k-CkTC-ki){</'k'-c_k'iCk'|) 



Z = / V(j)^V(j)i/Dck^Vck„e~ 



(2.3) 



where 



S = j dr 



where we defined 



( CkT c_kj, 

k 



'dr + e^ A 
A 9,-ek, 



CkT \ IAP 



(2.4) 



k 

From the action ( p.4|) , we can read off that in the mean-field approximation, the total free energy of 
the system is 

lAp 



^(|A|) = E(-^(^k) + -^(-i?k)) , 



(2.5) 



Since there is no universal convention as to whether the attractive interaction term should have a plus sign with a negative 
coupling constant V or a minus sign with positive V, we use \V\ which is unambiguously positive. 
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upon minimization with respect to |A|. Here, 

= -rin(i + e-^/^), (2.6) 

and 



We see the instabihty most clearly at T = 0, where -F(|A|) = -B(|A|). Then 

T{E) = e{-E)E, 



so 



I A 1 2 

E{\A\) - E{0) = N{0) I de{-^e^ + \A\^ " + W' ^^"^^ 

where A^(0) is the density of states at the Fermi level, and lod is the Debye frequency. Direct 
calculation shows that the integral behaves as |Apln|^ for |A| 0, whose non-analytic decrease 
will win over the analytic increase of the second term for small enough A, no matter how weak the 
attractive interaction \V\ is. By the same calculation, we can also see that the integral becomes 
analytic if we do not integrate e all the way up to zero, but to a finite negative energy. This means 
that the states close to the Fermi energy drive the BCS transition — they benefit most from opening 
of the gap |A|. Similarly, we shall see that the states at zero energy, that is the ZBSs, will drive the 
d ^ d + is transition. 

So far, the argument has shown the BCS instabihty only at T = 0. At finite temperatures, 

J^(E^) + J^(-E^) = -rin(2 + 2 cosh 
which, upon expansion in powers of |Ap, gives 

F(| A|) - F(0) = I A|2 - N{Q) I ^ tanh ^] + 0{\A\% (2.8) 
This shows that the system is unstable to the BCS transition at temperatures below T^. that satisfies 

1 _A.(0)/|tanh^ = 0. (2.9) 



In a similar way, we shall see below that T^, the transition temperature into the d+is state, is finite. 



III. THE D + IS INSTABILITY 



A. Formalism 



We now need to develop the formalism that will enable us to extend the strategy from Section 
H to the d + is case. We shall consider a single (2- dimensional) CuO plane, and model it by the 
Hamiltonian 

H = fd^rJ2 ^i(r)e(-«V)^,(r) + f rfVrfW(r - r')V^{(r)V^|(r')V^l(r')V^T(r), (3.1) 

where e is the band energy and V is the short-range interaction responsible for pairing. What makes 
this difficult problem tractable is the separation of energy scales (the Fermi energy Ep is much bigger 
than the superconducting gap A), which gives rise to separation of length scales Xp (Fermi wave 
length) and ^ (the coherence length). We may, therefore, expand in powers of the small parameter 
Xf/^', keeping the lowest non-trivial order is called the quasiclassical approximation. This procedure 
is usually done at the level of Green's function ]T^JTB| , which are thus transformed into Eilenberger 
functions that satisfy transport-like equations. 

Since we want to understand the d d + is transition in terms of quasiparticle eigenstates rather 
than Green's functions, we will perform this separation of scales at the operator level instead. We 
denote as 2A the width of the shell around the Fermi surface containing the states that take part in 
the pairing (see Fig. 2). We then factor out the fast Fermi-surface oscillations and define the slowly 
varying field operator 4'a,e{'r) fl^ by 

d'^k 



A 




^.,e(r)e*''^(^>^ (3.2) 
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FIG. 2. Fermi-surface decomposition of the Fourier transform. 



When we substitute this into (|0|), we obtain 



H 



E 



+ 



F.S. 



277 
" F.S. 

dkpie) dkpie') 

27r 27r 



(3.3) 



where "Vp{6) is the Fermi velocity at the point \<.f{6), and 



The derivation of (p.3|) is given in Appendix A. Note the hnearized kinetic energy in (|3.3|) , which 
will be crucial in the following. 

The Hamiltonian ( p.3|) gives rise to a partition function, which we can write as a path integral 
over the fermion fields ipafiiT^)- We can again decompose the interaction by the HS transformation, 
i.e., we can multiply the partition function by the constant 



dkpie) dkpie') 

27r 27" 



F.S. 

x(0e(r) -^teW^i~eW)(0e'(r) -^i-e'(r)^Te'(r))]- 



(3.4) 



In the mean-field approximation, the total free energy of the system equals the free energy given by 
the (single-particle) Hamiltonian 

V^Te(r) 



H = / rfV 



dk,{9).... ^ /vH^).(-.V) A,(r 



F.S. 
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A*(r) v^(^)-(zV)y V^i_,(r), 



F.S. 

upon minimization with respect to where we defined 



(3.5) 



A.(r)= / ^^V{e,e')Mr). (3.6) 

F.S. 

We shall write explicit formulae for the total energy and free energy in Sections ^ and |V| (formulae 
( [4.1|) and ( |5.1|) ). Here we just note that to calculate the single-particle contribution to the free 
energy, we will have to find the spectra of the Andreev Hamiltonians labeled by 6, i.e., we will need 
the energies Ee^n that satisfy |15| 



\^{9)■{-^V) A,(r) \//e,„(r)\ //.,„(r) 

= Ee,n . (3.7) 

A^(r) v^(^)-(^V)y \9eAr)J \9e,n{^) 

We note that the linear kinetic energy in ( p.3| ) makes this equation effectively one-dimensional, i.e., 
an independent equation for each line in the direction ^^{6). In the presence of the specularly 
reflecting boundary, we must find the Andreev spectra along refiected lines such as the one in Fig. 
la. Equivalently, we solve the equation on a straight line with the pairing potential A shown in Fig. 
lb. This is intuitively obvious; a derivation is given in Appendix B. 

As we mentioned in the Introduction, the spectrum along each trajectory having opposite signs of A 
at the two asymptotic ends will contain a zero-bias state. Its wave function is, up to a normalization 
constant. 



]exp{T J dp' A{e,p')), (3.8) 



.9{0,p)/zBS 

where the upper (lower) sign corresponds to A{9,p = — oo) < 0, A{9,p = +oo) > {A{9,p = 
— oo) > 0, A{9,p = +oo) < 0), so that the wave function is normalizable. In our notation, we 
will freely interchange the dependence on r (actually only on x, since the system is translationally 
invariant in the y-direction) with the dependence on the angle 9 and the coordinate p along the 
trajectory. Their relationship is obvious from Fig. 1. 
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B. Argument 



We now have all the tools needed to demonstrate the d + is transition in a way that brings out 
its physical mechanism. We follow the same line of thought as in Section ||: We go to the zero 
temperature, and look at the energy gains and losses when the s-wave component of A appears. 

For any s-wave pairing to appear, it is necessary that the part of the functional integral (|3.4|) over 
the s-component of converge, i.e., that V on top of the dominant ci-wave attraction contain also 
an s-wave part 0, 

v{e,9') = Vdie,9')-\Vs\. 

In this Subsection, we will show that this condition is also sufficient: At zero temperature, the system 
will favor the d + is state for an arbitrarily weak attraction Vg. 
With both d- and s-wave pairing present, 

00 (r) = (f)de{r) + (f)s{r). 

(The s-components of both V and are angle-independent.) We begin with the second term in 
( ^31) , which then is 

J Ztt Ztt 

F.S. 

J ZTT ZTT 

+ I>ll7 (3.9) 

F.S. 

where we used the orthogonality of the s and d components: 




(pde(r) = 0. 



F.S. F.S. F.S. 



We can also split up ( p.6|) into components and define 



^We use again |Vs| rather than 14. 
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F.S. 



A.(r) = / ^(-|K|)0.(r). (3.10) 

F.S. 

The (i-component of A was established well above Tg, so the change of the second term in ( 3.5|) due 
to the opening of a (small) s-wave gap will be 

'^■'-■"^ (3.11) 



just as in the BCS case. Due to the translational invariance in the y-direction, we will from now on 
write As(r) = As(x). Along the quasiclassical trajectory, x depends on both p and 9 (see Fig.l), so 
we will then write As{9,p). 

To examine the effect of the small s wave component on the quasiparticle energies, we need to 
look at the change of the spectra of the ID Andreev problems 

-tvpi9)d, A,i9,p)\ /Ui9,p)\ ffniO,p)\ 

]= Ee,n . (3.12) 

A,{9,p) zvF{9)dJ \9n{0,p)J \9n{e,p)J 

upon Ad{9,p) Ad{9,p) + As{9,p). As A^ is small, it can be treated as a perturbation; then the 

change of the quasi-particle energies to the lowest order is 

A,{9,p)\ (U9,pY 



Efl[As\= J dp{f:{9,p) g:{9,p)) 



-oo 

+ 00 



A:i9,p) / \9ni0,p) 
dpime, p)9ni9, p)Asi9, p) + g:i9, p)U9, p)A:i9, p)]. (3.13) 

Let us first look at the change of energy of the zero-energy bound states. Then from ( p.8|) 

9zBs{0,p) = TtfzBs{0,p), (3.14) 

+ 00 

4%s[^s] = ± / dp\f{9,p)\'2ImAs{9,p), (3.15) 

— oo 

where the upper (lower) sign corresponds to the +y- (— |/-)moving trajectory. We notice several 
things by looking at ( |3.15|) : 



so 
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• It depends only on ImAs, since ReAg just changes the position of the node in the total A{6,p), 
in which case the bound state remains at zero energy. Hence, we will assume ReAg = 0, and write 
A,{9,p)=zs{9,p). 

• Eg^zBsl^s] is non-zero due to the form of the bound-state wave function ( p.l4|) and due to the fact 
that s{6,p) does not change sign along the quasiclassical trajectory (by virtue of the s-symmetry). 
Out of the two possibilities for the sign of s, we will choose s{6, p) > in the following, which means 
all the -|-y-moving states are shifted up in energy, whereas the — y-moving states are shifted down. 

Since we are at zero temperature, only the states that move down from zero energy will be occupied. 
We can then argue similarly as in the BCS case: opening of the additional s-wave gap costs the 
system energy s^/|Vs| (from ( 3.11 )) but the quasiparticles save energy ~ s. The lowering of the 
quasiparticle energy is only linear in s, i.e., not as dramatic as the non-analytic decrease in the BCS 
case, but nevertheless it beats the quadratic increase for small enough s. Thus, for an arbitrarily 
small but non-zero interaction IKjI, s = cannot be a minimum of the total energy, and the additional 
s-wave gap phase shifted by 7r/2 from the rf-wave gap will appear. From the formula (|3.15|) , we see 
that the superconductor will benefit from opening up the gap only close to the surface where |/p 
is effectively non-zero, so the transition into the d + is state is a surface effect. The decay into the 
bulk will be discussed more quantitatively in the next section. 

We should also note that the remaining states on the quasiclassical trajectories do not change this 
situation, that is, they do not contribute linearly to the change of the total quasiparticle energy. 
Due to the time-reversal symmetry in the pure (i-wave state, every state on a given quasiclassical 
trajectory corresponds to a state of the same energy on a reversed trajectory. Indeed, if we label the 
coordinate along the trajectory reversed to the one in (|3.12| ) as p = — p, then the Hamiltonian on 
the reversed trajectory is 

-tvF{-d)d/, A,{-e,p)\ f ivF{e)d, A,(e,-p) 
Ad{-erp) ivF{-~9)d~J \Ad{e,-p) -ivF{e)dp 

since vf{—0) = vf{0), so 

fn{-e,p)\ _ /9n{e,-p) 

.9n{-9,p)) [fnie,~p) 
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will also have energy Eg^n- Now (|3.13| ) implies that to the first order, a small ImAs will shift the 
energies of the two corresponding states by an equal amount with opposite signs. Hence, the only 
way they can linearly contribute to the total energy at T = is when one of them crosses zero and 
thus changes its occupancy, which happens only when their original energy (in absolute value) is 
smaller than the s-wave gap. But as s — 0, there will be fewer and fewer such states in smaller 
and smaller neighborhoods of the c/-wave nodes. It is only the ZBSs that change their occupancy for 
arbitrarily small s. We thus conclude that the onset of the transition into the d + is state is driven 
by these states. 

IV. S-WAVE GAP AT r = 

In the study of the instability of the pure d state in the last section, we used first-order perturbation 
theory since s ^ at the onset of d+is. Now we will argue that this theory holds up to the the actual 
value of s, i.e. s << \Ad\. As we will show in Section 0, s(T = 0) ~ Tg, the transition temperature 
into the d + is state. Because ~ 7 K, it is much smaller than T^, the superconducting transition 
temperature of the order of lOOK, which sets the scale for A^^. Figure 3 shows the magnitude of the 
two gaps as a function of angle around a quarter of the Fermi surface. We see that the required 
inequality s << \Ad\ holds for most of the Fermi surface except for small neighborhoods of the 
nodes. First-order perturbation theory certainly breaks down there, but upon averaging over the 
Fermi surface, the nodes will only introduce an error of the order Ts/Td. Thus, we will use that theory 
to obtain s{x) at T = 0. As discussed at the end of Subsection [IIIB|, first-order perturbation theory 



implies that we have to look only at the zero-energy states. Also, since s is a small perturbation, we 
shall neglect its effect on A^^. 
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^0 



s(T=0) 



e 

-Ji/4 Ji/4 

FIG. 3. The magnitude of the d- and s-wave order parameters around the Fermi surface. 

Now we can write down the energy due to s per unit length of the surface (the y-direction) as a 
functional of six): 



9e(-7r/2,0) 



27r 



Ee[s{d,p)] cos 6, 



(4.1) 



where Eg[s\ is given by ( 3.15 ); for the rest of this section, we shall drop the superscript "(1)", since 
we shall be using only the first-order formula. We freely interchange s{x) for s{6,p); the relation 



between the two is discussed below ( p.ll|) . Note the correct dimensions: the x- integration makes 
from energy per unit area into energy per unit length. In the second term, the integrand 
is energy and the dimension of the measure is kp, i.e., inverse length. The extra factor of cos^ in 
the second integral accounts for the difference of the density of trajectories along the y-direction 
compared to their angle-independent intrinsic density (measured perpendicularly to their direction), 
as shown in Fig. 4. In the second term in ([4.1| ), we sum up only the occupied —?/- moving states for 
which Eq < according to ( p.l5|) . 
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FIG. 4. The decrease of the density of the trajectories in the y-direction by the factor cos6'. 

We obtain s{x) by minimizing ([4.1|) . Let us first make an order-of-magnitude estimate 

E[s] ~ - kps, (4.2) 

since s will extend into the bulk only as far as the coherence length ^ = hvp/Ao (A^ is the amplitude 
of the d wave), and from ( p.l5|) , we see Eg[s] ~ s. The angular averaging will, up to numerical factors 



of order unity, multiply Eo[s] by kp. Minimization of ( |4.2| ) will give 

kplVsl 



(4.3) 



By taking s ~ ImeV from the experiment, kp ~ 1A~^, and ^ ~ lOA, we get an estimate for the 
strength of the s-wave pairing 

IKI ~ lOmeVA^ 



We minimize ( [4.1| ) exactly by solving 

5E[s] 



Ss(x) 



ee(-7r/2,0) -oo 



Now 
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5s(6,p) ^, ^, ^, 

o[x — p cos u) + o{x + p cos ( 



^(p-^)+^(P+^))^ (4-4) 



5s{x) 

cos 6 V"^'" cos 6'^ ' "\r ' (^Qg^' 

since cos^ > 0, and p, unlike x, can be both positive and negative. The factors of cos^ cancel, and 
we obtain 

s{x) = m I ^(iM^)r+im— ^)r). (4.5) 

J Ztc \ cosd COS 6 / 

6le(-7r/2,0) 

Physically, we get two terms on the right-hand side because for each angle 6, there are two trajectories 
contributing to s at a given point as shown in Fig. 5. 



FIG. 5. Two contributions from the same 9 to the pairing potential at the point x. 



We should remark here that we also obtain the formula (|4.5| ) when we calculate the contribution 
from the occupied (— y-moving) bound states to the pairing potential in the gap equation. This is 
done in Appendix C. The result is 

A.M,,,.,iv;i / ^J!i^nne,-^)f^\m-^m (4.6) 

9e{-f,0) 

in agreement with ( [4.5| ). This formula, however, shows more clearly the internal consistency of the 
picture: For in Fig. 1, the additional is potential pushes down the — y-moving states if s > 0. 



As ( [4.6|) shows, these states, in turn, give rise to = i x positive. 

We should note here that A^ is absent on the right-hand side of (|4.6|) , so the gap equation in this 
case (unlike in the BCS theory) is an explicit formula for the gap. The physical reason for this is 
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that s{x) is considered small, so we neglect the change of the bound-state wave functions due to its 
presence. The only effect of s{x) we are taking into account is the change of the occupancy of the 
zero-energy states, which, by (|3.15| ), depends only on the sign of s, not on its detailed shape. This 
is why s{x) does not feed back into the right-hand side of (|4.6| ). 

To estimate the decay of into the bulk, we shall assume to be constant in space and with 
the angular dependence 

Aa,e{r) = AoSm29, (4.7) 

which should hold for 

hVF 

Also, we shall assume a spherical (circular) Fermi surface, 

dkpie) = kpde. 

Then the wave function of the — y-moving bound states, including the normalization, will be 

'f(^^p)\ I\sm20\ /I 



-|psin2e 



so 



kplVsl } de \sm2e\ .^i.^^^^el 

six) = / — 2 X e €icosfl I = 

^ ' f J 2n 2 

-tt/2 

^-^l^sl f rn ■ n n -4sinef r,\ 

dOsmOcosOe « . (4.9) 







We can do the integral by substitution sin 6 = t, which gives 



2\ n 1 

4*2 



(4.10) 



We can neglect the contribution from the upper limit because it is effectively non-zero only for 
X < ^/4, where our assumption of constant A^ does not hold. The lower limit should have been 
at Ts/Td, rather than at 0, to exclude the trajectories close to the nodes where the first-order 
perturbation theory breaks down. That cuts off the lower-bound contribution at a; ~ ~ lOOA, 
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beyond which we would need a more refined theory for the behavior of the quasiparticles around the 
nodes. For x much smaller than this distance, we can neglect the first term on the right hand side 



of ( [4.10|) , and replace the exponential by 1. We conclude, therefore, that 



for 



i<x« — e 



We see that s = kp\Vs\/^ times a function that is of order unity for x < ^, and decays fast for x > ^, 
as expected. 



V. TRANSITION TEMPERATURE 

So far, we have shown the instability d ^ d + is only at T = 0. Just as in the BCS case, it remains 
to be demonstrated that the transition temperature Tg is finite. We therefore must study the free 
energy of the system, which we obtain from (|4.1j ) when we replace Eg[s] by J-'{Eg[s]), the free energy 
of a single fermion level (see (|2.6|) ), that is, 

F[s] = J dx'^ + kF I ^cos^(-T)ln(l + e-^«H/^). (5.1) 

J \VA J ZTT 

' ' -7r/2 

Minimization of this functional will give an equation for s{x) that again agrees with the contribution 
to the gap equation from the ZBSs. As we see from ( ^.1| ), the variational equation for s will now 
be very non-linear; it will no longer be an explicit formula for s. The reason is that at finite 
temperatures, the occupancy of a given state depends on the value of its energy. Even in first-order 
perturbation theory, this value depends on the shape of s{x), not just its sign, so s{x) enters through 
the Fermi function into the right-hand side of the gap equation, making it non-linear and therefore 
difficult to solve. 
We still can make an order-of-magnitude estimate of F as follows 

ln(l + e-^«W/^) + ln(l + e'^-^W/^) = ln[(l + e"^«H/^)(l + e^^'W/^)] 
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ln(2 + 2cosh.^^f'] 



T ' 

-1114 + — + 0(5'), 



SO 



i^M-i^[0]~^^(^-^)+O(s^). (5.2) 

From ( |5.2| ) we see that the system is unstable to the transition to the d ^ is state below the 
temperature ~ kp\ys\li, which is therefore of the same order of magnitude as |As|r=o. (see 

D). 

Following we can trade the coupling constant K for the transition temperature, T^s of a BCS 



superconductor with this coupling, T^s ~ e ^/I^^L Then 

-1 



InTcs 

Hence, Tg increases sharply close to T^s = 0, which is consistent with the numerical results P,p!0 



(5.3) 



VI. CURRENT 

To study the surface current in the d+is state, we shall go back to T = for simplicity. We observe 
that the states on the +?/-moving quasiclassical trajectory from Fig. 1 will, upon the transition into 
the d + is state with s > 0, feel the pairing potential shown in Fig. 6. As p goes from —00 to 
+00, the twist of the phase of the order parameter is clockwise (from vr to 0) for an +?/-moving 
trajectory and counterclockwise (from to tt) for a — ?/-moving one. In both cases, this implies 
current flowing in the —y direction. This agrees with our previous calculations that showed that the 
—y- (+y-)moving bound states will be (un)occupied if s > 0. 
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ImA 



Re A 



FIG. 6. The pairing potential along the trajectory in Fig. la). The corresponding twist of the phase of the order parameter 
is clockwise. 

The agreement is quantitative, as we can easily check. In the hnearized Andreev formahsm, the 
contribution to the current density from a given state is 

£''\0,p)=evMnie,p)\' + \g^ie,p)\% (6.1) 

that is, the charge of the state times its (Fermi) velocity times the occupation of that state. In 
our case, all of the current is carried by the occupied bound states because the contributions from 
the remaining pairs of corresponding countermoving states cancel each other out (see the end of 
Subsection III B ). To calculate the total current density in the ?/-direction, we again have to include 



the contribution from both the incoming and the outgoing part of each — y-moving trajectory (see 
Fig. 5), and we have to project onto the y-direction 



-tt/2 



(6.2) 



On the other hand, in terms of the one-dimensional density n^^'^'^ = kp/Tr and the phase of the 
order parameter ip{0,p) along the trajectory, 

= ^-^dMO,p), (6.3) 

since for a spherical Fermi surface vp = kp/m. The formula for the total surface-current density will 
be the same as (|6.2|), except that we now have to integrate over both +y- and — ^/-moving trajectories 
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-7r/2 



COS^ 



(6.4) 



We do not expect the current densities ( |6.2| ) and ( |6.4| ) to be the same at a given point because the 
formula ( |6.3|) has corrections, which are higher-order derivatives of (p. Those corrections will not, 
however, contribute to the total surface current 

ly = dxjy{x), (6.5) 

which should then come out the same in the two calculations. Indeed, the bound states give us 

+0O 

{IzBs)y = kF j —sin 6* COS 6* j dpj^zBsi^y p)^ 



(6.6) 



smce 

+00 



due to the normalization of the wave functions; the minus sign indicates that the current is flowing 
in the —y direction. The formula for the total current in terms of the order-parameter phase is 

iIo.p.)y = kF I — sin^cos^ I dpi'^^ie,p). (6.7) 



-7r/2 



Now 



dp£^{e,p) = ^(^(^,+oo) -^(^,-oo)) 
evp 



TT 

oo 



so 



sgn(^), (6.8) 



iIo.p.)y = = {IzBs)y, (6-9) 



since sgn(6') sin 6' is an even function, so the factor 1/2 in ( |6.8|) compensates for the doubling of the 
integration domain of 6 in (|6.4| ) compared to ( |6.2|) . To get an order-of-magnitude estimate, we put 
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e ~ 
Vp ~ lO^m/s 

and get \Iy\ ~ 10~^A per CuO plane. From the approximate form of the bound state wave functions 
introduced in the previous section, we can also estimate the spatial distribution of the current density 



jy{x) = Aevphp I — sin 6' 
J 2% 



f 



X 



-Tr/2 
tt/2 

AevpkF f dd 



cos 9 



2evpkp 



J 27r 

dtfe~T\ (6.10) 







By the same argument as presented in the last section, we find that for ^ < x « 

The extra power of x in the denominator in ( |S.11D , compared to ( [4.11 ), comes from the directional 
sine in (|6.2|). 



The surface current induces magnetic field, which will be screened by the diamagnetic current in 
the superconductor. The total current density therefore is 

{jtot{x))y = {jzBs{x))y + {jdm{x))y (6.12) 

According to ( |S.11| ), the current is localized within the distance ~ ^ from the surface, which is much 
smaller than A, the in-plane penetration depth, because the cuprates are strongly type-2 supercon- 
ductors. Hence, the diamagnetic response does not resolve the internal structure of {jzBs{x))y, and 
we can estimate 

{jdm{x))y = (jd,n(0)),e-"/\ (6.13) 
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FIG. 7. Side view of the ah planes with the current flowing in the —y direction. The induced magnetic field is along the 
c-axis. 

The surface current will be screened completely, because the magnetic field it induces is smaller 
than B^^, the lower critical field in the c-direction. Indeed, even if the current fiows in the same 
direction along all the CuO planes (as shown in Fig. 7), the magnetic field at distances x > dc (the 
interplane spacing ~ lOA) from the surface will be 

B = ^1 ~ lO^G, (6.14) 

c dc 

which is smaller by an order of magnitude than B^^ for YBCO (see |T^). The complete screening 
implies 



{Itot)y = I dx{jtot{x))y = 0, (6.15) 

which together with (|6.9|) gives 







iUmmy = (6.16) 



For ^ < X << jT^, we can approximate the exponential in ( |6.13| ) by 1, so 
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{jtot{x))y 



(6.17) 



167r 

This changes sign at distance 

Xo ~ ^k'/\ (6.18) 

where k = X/C, is the Ginzburg-Landau parameter. Due to the one-third power, Xo ~ ^ for reasonable 
values of k (say, between 50 and 500). This is consistent with the numerical results 0. 
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Note that the two calculations of the surface current agree (see (|6.9| )) because the ZBSs moving 
in the direction of the current are shifted down in energy and thus occupied, whereas those moving 
against the current are shifted up and unoccupied. We wish to stress that this is exactly opposite 
to the sign of the Doppler shift: the states moving along the current would be Doppler-shifted up, 
whereas those moving against the current would be Doppler-shifted down. 

This point is further supported by the analogy between the ZBSs and low-lying excitations in a 
core of an s-wave vortex. We consider an idealized case: A = inside the vortex (at distances from 
the center smaller than R), and |A| = const outside with the phase winding counterclockwise once 
around. We look at a quasiclassical trajectory passing close to the center of the vortex. We denote 
the coordinate along the trajectory as p again and the phase at the intersection point with the vortex 
edge as ip±, see Fig. 8. Then the energy of a low-lying excitation moving from p = — oo to p = -|-oo 
on that trajectory is jTsI 



For a trajectory passing through the center, y9+ — y9_ = vr, so = 0, and the low-lying excitation 
is a ZBS. If we now shift the trajectory slightly to the left as shown in Fig. 8, then ip^ — tt 
and E > 0. In a real vortex, A would be non-zero even inside, and the phase of A would wind 
clockwise as we go from p = — oo to p = +oo, so we are going against the current. Moreover, 
(p(p = +oo) — ip{p = — oo) = TT mod 2tt, so A behaves the same way as in Fig. 6. Hence, the 
d ^ d + is transition is analogous to shifting the quasiclassical trajectory away from the vortex 
center. In both cases, the ZBS will have a positive energy if it is moving against the current and 
negative energy if it is moving in the direction of the current. 




(6.19) 
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FIG. 8. A quasiclassical trajectory going through the core of an s-wave vortex. 

VII. DISCUSSION 

In order to understand the basic physical mechanism of the d ^ d + is transition, we considered 
the change of the free energy of a CuO half-plane when an s-wave component of A appears close to 
the 110 surface. By the Hubbard- Stratonovich transformation and the mean-field approximation, 
we decomposed the total free energy into the contribution from the single-particle states, which is 
decreased by Im A^, and the HS term, which increases quadratically with A^. Using first-order 
perturbation theory, we saw that the system favors the d + is state at T = 0, and that the transition 
is driven by the zero-bias states. Based on this argument and on the separation of energy scales 
associated with the d- and s-wave components of A, we then calculated A,, at zero temperature and 
estimated the transition temperature. Finally, we discussed the surface current in the d + is state. 
We saw that it is carried by the occupied ZBSs; these states are not Doppler-shifted by the current. 
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APPENDIX A: DERIVATION OF THE QUASICLASSICAL HAMILTONIAN 

We derive ( |3.3|) from (|3.1| ) by expressing the original field V'<T(r) in terms of the slowly varying field 
ip^g^r), see {\i.'2\). We first substitute into the kinetic-energy operator: 



e(-iV)^.(r) = / ^^e^^-W-^e(k^(^) - zV)^.,,(r) = 

F.S. 

J dk^^^^^ey. ^.i^^^^^e)) + Vke(k)|k,(^^) • (-^V) + ■ ■ ■] V^.,,(r) ^ 

'■v^(0)-(-2V)^.,,(r), (Al) 



F.S. 



F.S. 

since 



e^kp) = 0, andVke(k)|k^ = vp. 

The kinetic energy therefore is 

J 1, I J „ J 2tx J 2tt 



F.S. F.S. 



xe.'(r)v^(^)-(-2V)^.,,(r). (A2) 



Now e*'^''^'^^)^''^^^'^)''' oscillates with a wavelength much shorter than the lengthscale on which ^p^^g(r) 
changes. Thus, the integral will be zero unless ki?(^^) — \<.f{0') = 0, so we can effectively drop one 
integration over the Fermi surface, and obtain the kinetic energy of the form 

/^'-E / ^4A-)M0) ■ HV)^.Ar). (A3) 

The potential energy is now given by 

,2 ,2 >dkF{ei) dkF{92) dhpiOs) dkpie^) 



d rd r Vir — r x 

27r 27r 27r 27r ^ ' 



X 



exp[^(-k^(^^i) • r - k^(^^2) ■ r' + \^f{9,) ■ r' + kF{e,) ■ r)] x ^\,^{r)^pl^{r')^P^g,{r')^^g,{r 
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~ f d2^^2^, dkF{0l) dkpm dkpm dkF{9i) ^ ^n^.(k^(e,)-kWe3)) (r-r') ^ 

J 2-K 2-K 2-K 2-K ^ ' 

xWe,(r)We.(r)^ie,(r)^Te4(r) x exp[z(-k^(ei) - 1.^(^2) + k^(^3) + k^(^4)) • r] 
f 2 ,2 ,dkF{9i) dkF{92) dkpiOs) dkpiO^) 

= 7^^^^^^^ ^m,^3)x 

xV^{,^(r)^{,^(r)^i,3(r)^T,,(r) x exp[z(-k^(0i) - k^(^^2) + k^(^3) + k^(^4)) • r], (A4) 

where we used the assumption that V changes on a much shorter length scale than '?/'o-,6»(r), performed 
the r'-integration, and introduced the Fourier transform 

Since V{t) = V{—r), we see that V{9,9') = V{9',9). Again, the integral vanishes unless the sum of 
the four momenta is zero. Out of the various ways that this may happen, we pick only the one that 
contributes to the singlet pairing, namely 

kp{9,) + kp{92) = kF{9s) + k,.(^4) = 0, t.e. 

9^ + 92 = 9-^ + 9^ = mod 2%. (A5) 

Since we now have two constraints, we can drop two Fermi-surface integrations, and obtain the 
potential energy of the form 

J Zn Zn I + 

The total Hamiltonian in the quasiclassical approximation then is 

H = / ^<,(r)vH^) ■ (- W.,.(r) + 

°" F.S. 



(A6) 

Z7I Z7I ' ■ J 

F.S. 

The derivation of ( |A6| ) from ( ^.1| ) is far from rigorous. We could give a somewhat better, although 
much longer, argument. We believe the approximations used here are equivalent to the approxima- 
tion in the Eilenberger formalism, because the Eilenberger equations can now be rigorously (apart 
from the mean- field approximation) derived from (|A6|) . 
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APPENDIX B: BOUNDARY CONDITIONS AT THE SURFACE 

We show the effect of the boundary on the Andreev spectrum. In general, the boundary will 
cause mixing of different 6''s. For each 6', though, we have a different Andreev equation, so adding 
the solutions of ( |3.7|) for different 6''s does not make sense. However, the Andreev wave functions 
describe only the slow variation of our excitations (changes on the length-scale ^). The full wave 
functions containing the rapid oscillations as well are 

e^k^W-r^ (Bl) 

and these describe the single-particle excitations of the same Hamiltonian (|3.1|) (in the mean-field 
approximation), so those can be added. If we assume a specularly reflecting boundary, then the 
wave function will contain only two terms: 

>^e,„,n(r) 



such that 

Oin + Oout = TT mod 27r, (B2) 

since the angles are measured from the positive-x semiaxis, see Fig. la. The Dirichlet boundary 
condition gives 

f0out,n{r) \ ( fei„,n{r) 

resurface ' 



whereas the Neumann boundary condition gives 
since 



resurface 



n-{kFi9in)-kp{9out))=0 



for n perpendicular to the surface, and we used 



7\ //' 

kp ^ \ » n Vr \ 
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to neglect the gradient of the Andreev wave function. Since (|3.71) is hnear, we can drop multiphcative 
constants and simply assume 

/0ont,n(r) \ / /e„,n(r) \ 

(B3) 

9o„t,n(r)/ \fl'6i,„,n(r) / 

at the surface for either choice of the boundary condition. As 6'j„ is uniquely determined by Oout 
through the relation (|B2|) , we shall label the potential A along the trajectory as well as the solutions 
of the corresponding Andreev equation by 9 out- We shall drop the subscript "owt" everywhere except 
in Appendix C, where we will need to distinguish 9 out) the label for a trajectory as in Fig. 1, from 
9, the label for a position on the Fermi surface as in Fig. 2. 

APPENDIX C: GAP EQUATION 



We obtain the gap equation by substituting for (t)e{j) in ( p.6|) its mean-field value, that is, the pair- 
ing amplitude ^^(r) = (?/'|_e(r)-?/'|0(r)). To calculate this amplitude, we expand the field operators 
into energy eigenstates 

^T^(r)\ ^ (U9,p)\ 

= E7.,n . (CI) 

.^-.(r); n \gn{9,p)) 

Equation ( Pl| ) gives at T = 

{^P^.eir)^P^oir)) = UO, p)9*A9, p){^l^ne,n) = 

n,n' 

= J2ei-Ee,n)fn{9,p)g:{9,p). (C2) 



Note that in we explicitly sum over both positive and negative energies, unlike the Bogoliubov- 
de Gennes (BdG) formalism where we can sum over positive energies only, using the fact that 



and (C3) 
t;„(r); \Kiv)l 

are both solutions of the BdG equations with energies equal in absolute value and opposite in sign. 
However, this symmetry is lost here because the Andreev wave functions corresponding to the BdG 
wave functions ( |C3|) live on different quasiclassical trajectories. 
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Close to the surface, we have to remember again that each hne contributes to the pairing amphtude 
at a given point for two directions 6*; see Fig. 5. We will, therefore, have to distinguish between the 
label of the trajectory Oout and the label for the pairing amplitude 6. Specifically, 

TT 

9out = for 6'e(--, 0) and 

eout = ~7r-0 for ^e(-7r,-|), (C4) 
so the contribution to the pairing amplitude from the — ^/-moving bound states will be 

= fiO, ^)9*iO, for 0e(-^, 0) (C5) 

cos 6 cos 6 2 



and 



{'ipi-e{j)%l)^e{j))zBS = f{Oout, ^ — )9*{6out, ^ 



cos 9out COS 9out 



Substituting ( PSP and ( |C6D into ( p.6| ) gives 



e'e(-n,0) 



J ZTT COS t7o„j 



+ V{9, -TT - (C7) 

COS dout 



where 

mcj = -iKi+w,o, 



and we used ( p.l4 ). The contribution to the s-wave component of the pairing potential from the 



occupied bound states therefore is 

^s{x)zBS = r\Vs\ f + (C8) 

J 271 cosO'f cosO'i- 

For the calculation of the d-wave component of Azbs, we will assume that the unperturbed is 
antisymmetric around its vertical node, 
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Ad,,(r) = -Arf,_^_,(r). (C9) 
Presumably, A^ arises from an antisymmetric interaction 

Vd{e,e') = -Vd{e,-T:-e'). (cio) 

Along the quasiclassical trajectory, ([C9| ) means 

A,(^^,p) = -A,(^^,-p), (Cll) 

which, by ( ^78|) , imphes 

IMP)P = IM-P)r. (C12) 

Thus, under these assumptions, 

^,,e{x)zBS = {-^) I '^^^^me, + V,{e, -vr - ^^)r = 

= 0. (C13) 
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